SoSe2022

Folienübersicht

Poweranalyse und Berechnung des Stichprobenumfangs

Data Science 2 | Folienübersicht

Ãœbliche Experimentfragen

Nachdem das Design des Experiments und die statistische Analyse bestimmt wurden, ergeben sich oft folgende 3 Fragen:

  1. Welcher Stichprobenumfang ist nötig?
  2. Wie hoch ist die Teststärke oder Power des Experiments?
  3. Was ist der kleinste Effekt, der mit dem Experiment nachgewiesen werden kann?
Data Science 2 | Folienübersicht

Einflussgrößen

Für die Beantwortung dieser Fragen braucht es folgende 5 Informationen:

  1. Stichprobenumfang (N)
  2. Effektgröße (ES oder \(\delta\)) → die Mindestgröße, die nachweisbar ist oder die vorhergesagte Größe des Effekts.
    • Je nach Experiment handelt es sich um die Differenz zwischen 2 Mittelwerten (t-Test), Unterschiede in Proportionen, das Verhältnis zwischen ‘between-group’ und ‘within-group’ Variabilität (ANOVA) oder die Korrelation zwischen 2 Variablen.
  3. Variabilität des zu messenden Merkmals (\(\sigma^2\)) → bezieht sich beim t-Test auf die gepoolte Standardabweichung, bei der ANOVA auf die ‘within-group’ Variabilität.
  4. Teststärke (Power) → wenn diese fix sein soll, wird oft Wert von 0.8 oder 0.9 gewählt.
  5. Signifikanzniveau (\(\alpha\)) → gewöhnlich festgelegt auf 0.05.
Data Science 2 | Folienübersicht

Zusammenhang

Die universale Power-Gleichung

\[Power \propto \frac{ES~\alpha~\sqrt{n}}{\sigma^2}\]

Die Teststärke steigt proportional* mit

  • steigender Effektgröße
  • steigendem \(\alpha\) (wenn das Signifikanzniveau weniger strikt wird)
  • zunehmendem Stichprobenumfang
  • abnehmender Varianz

→ Durch Umstellung der Gleichung kann jede dieser 5 Größen berechnet werden, wenn die anderen 4 gegeben sind.

*Wie genau die Teststärke sich ändert hängt von der Analyse ab.

Data Science 2 | Folienübersicht

Zum Glück gibt es auch R Funktionen

  • Je nach statistischer Analyse gibt es eine spezifische Funktion mit der die Stichprobengröße (n), Teststärke (power) und testbare Effektgröße (delta) ermittelt werden kann, z.B:
    • t-Test: power.t.test()
    • Zweistichprobentest für Proportionen: power.prop.test()
    • 1-faktorielle ANOVA: power.anova.test()

Beispiel t-Test 

power.t.test(n = NULL, delta = NULL, sd = 1, sig.level = 0.05,
             power = NULL,
             type = c("two.sample", "one.sample", "paired"),
             alternative = c("two.sided", "one.sided"),
             strict = FALSE, tol = .Machine$double.eps^0.25)
Data Science 2 | Folienübersicht

Berechnung des Stichprobenumfangs

t-Test

N beim unabhängigen t-Test

Um die Stichprobengröße für einen unabhängigen t-Test zu ermitteln, brauchen wir einfach nur die Formel der Teststatistik nach n umwandeln. Es gilt für gleiche Varianzen und Stichprobengrößen:

\[t=\frac{(\bar{X}_1-\bar{X}_2)}{\sqrt{2\frac{s_p^2}{n}}}=\frac{d}{\sqrt{2\frac{s_p^2}{n}}}~~\Rightarrow~~\sqrt{2\frac{s_p^2}{n}}=\frac{d}{t}~~\Rightarrow~~2\frac{s_p^2}{n}=\frac{d^2}{t^2}~~\Rightarrow~~n = \frac{2t^2s_p^2}{d^2}\]

  • \(s_p\) ist die gepoolte Standardabweichung bei gleicher Stichprobengröße: \(\sqrt{\frac{s_1^2+s_2^2}{2}}\)
  • \(t\) hängt von der Wahl der Power und des Signifikanzniveaus ab.
  • Werden diese auf 1-\(\beta\) = 0.8 und \(\alpha\) = 0.05 gesetzt, so beträgt t = 2.8 und \(t^2 = 7.84 \approx 8\) (die Quantilen einer Normalverteilung bei 0.05/2 und 0.8 werden addiert).
  • Somit gilt: \(n = \frac{2*8*s_p^2}{d^2}\)
Data Science 2 | Folienübersicht

Berechnung des Stichprobenumfangs

t-Test - Beispiel 1

Beispiel Zugverhalten: Buchfinken vs. Mönchsgrasmücken

Kenngröße Buchfink Mönchsgrasmücke
Mittelwert 1800km 3000km
Standardabweichung s ±900km ±1000km
Stichprobengröße n 20 30


Buchfink


Mönchsgrasmücke
d <- abs(1800-3000)
s_p <- sqrt((900^2 + 1000^2) / 2)
(n <- (16*s_p^2) / d^2)
[1] 10.05556
  • → Um einen Unterschied von 1200km statistisch nachweisen zu können, braucht es mindestens 10 Replikate pro Gruppe.

Bildquellen: Wikipedia - Buchfink (CC BY-SA 2.5) Mönchsgrasmücke (CC0)

Data Science 2 | Folienübersicht

Berechnung des Stichprobenumfangs

t-Test - Beispiel 1

Beispiel Zugverhalten: Buchfinken vs. Mönchsgrasmücken

Kenngröße Buchfink Mönchsgrasmücke
Mittelwert 1800km 3000km
Standardabweichung s ±900km ±1000km
Stichprobengröße n 20 30


Buchfink


Mönchsgrasmücke
d <- abs(1800-3000)
s_p <- sqrt((900^2 + 1000^2) / 2)
(n <- (16*s_p^2) / d^2)
[1] 10.05556
  • → Um einen Unterschied von 1200km statistisch nachweisen zu können, braucht es mindestens 10 Replikate pro Gruppe.
power.t.test(power=0.8, sd=s_p, delta=d)
     Two-sample t test power calculation 

              n = 10.9148
          delta = 1200
             sd = 951.3149
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

Bildquellen: Wikipedia - Buchfink (CC BY-SA 2.5) Mönchsgrasmücke (CC0)

Data Science 2 | Folienübersicht

Berechnung des Stichprobenumfangs

t-Test - Beispiel 2

Beispiel Zugverhalten: Buchfinken vs. Mönchsgrasmücken

Und wie groß müsste die Stichprobe sein, wenn wir bereits einen Unterschied von 600km bzw. 300km statistisch nachweisen wollen?

delta von 600km 

power.t.test(power = 0.8, sd = s_p, delta = 600)
     Two-sample t test power calculation 

              n = 40.44558
          delta = 600
             sd = 951.3149
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

delta von 300km 

power.t.test(power = 0.8, sd = s_p, delta = 300)
     Two-sample t test power calculation 

              n = 158.8158
          delta = 300
             sd = 951.3149
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group
Data Science 2 | Folienübersicht

Shiny App zur Berechnung von N

Unabhängiger t-Test

Data Science 2 | Folienübersicht

Stichprobenkalkulation vs. Poweranalyse

Tuning der Stichprobengröße

  • Faktoren, die den Stichprobenumfang bestimmen:
    • Durchführbarkeit / Zeit / Konvention
  • Dadurch ist eigentlich oft schon im Vorfeld der Umfang festgelegt.
  • Aufgrund der oft überraschend hohen Stichprobengröße wird dann gerne so lange an den andere Stellschrauben gedreht, bis das erwünschte N raus kommt → das ist aber nicht Sinn und Zweck der Analyse.
  • Die Stichprobenanalyse wird daher am besten genutzt rauszufinden, ob das geplante Experiment Sinn macht oder nicht:
    • Wenn die Teststärke beim geplanten N nur noch bei 40% liegt, ist es unwahrscheinlich ein signifikantes Ergebnis zu erhalten - auch wenn es einen Effekt tatsächlich gibt.
  • Anstatt die Power festzulegen und N zu bestimmen, kann auch N auf einen maximal durchführbaren Wert festgelegt werden, um dann die Power zu berechnen.
Data Science 2 | Folienübersicht

Poweranalyse

t-Test - Beispiel 1

Beispiel Zugverhalten: Buchfinken vs. Mönchsgrasmücken

N = 20 

power.t.test(n = 20, sd = s_p, delta = d)
     Two-sample t test power calculation 

              n = 20
          delta = 1200
             sd = 951.3149
      sig.level = 0.05
          power = 0.9729651
    alternative = two.sided

NOTE: n is number in *each* group

Im vorliegenden Experiment betrug der Stichprobenumfang 20 und 30. Angenommen in beiden Fällen war N = 20, wie hoch ist dann die Teststärke gewesen?

N = 5 

power.t.test(n = 5, sd = s_p, delta = d)
     Two-sample t test power calculation 

              n = 5
          delta = 1200
             sd = 951.3149
      sig.level = 0.05
          power = 0.4190796
    alternative = two.sided

NOTE: n is number in *each* group

Wie hoch wäre die Teststärke bei 5 Vogelmessungen gewesen, um einen Unterschied von 1200 km statistisch nachweisen zu können?

Data Science 2 | Folienübersicht

Shiny App zur Berechnung der Power

Unabhängiger t-Test

Data Science 2 | Folienübersicht

Berechnung des minimal nachweisbaren Effektes

t-Test - Beispiel

Angenommen es soll die Untersuchung des Zugverhaltens wiederholt werden, es können aber nur 8 Vögel jeweils untersucht werden. Welche Differenz könnte dann überhaupt statistisch nachgewiesen werden?

N = 5 

power.t.test(power = 0.8, n = 8, sd = s_p)
     Two-sample t test power calculation 

              n = 8
          delta = 1433.285
             sd = 951.3149
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

→ Die Differenz müsste mindestens 1433.3 km betragen, um überhaupt einen signifikanten Unterschied zwischen den beiden Arten nachweisen zu können.

Data Science 2 | Folienübersicht

Sog. ‘Power Curves’

1

  • Anstatt immer nur einzelne Szenarios durchzuspielen ist es hilfreich, die Teststärke für mehrere Kombinationen an Parametern als Kurve zu visualisieren:
# Erstelle df mit versch. Kombinationen
df <- expand.grid(n = seq(5, 40, 2), sd = c(500,750,1000, 1250), 
  delta = c(500,1000), power = NA)
dim(df)
[1] 144   4
# Powerberechnung für jede der Kombinationen
for (i in 1:nrow(df)) {
  df$power[i] <- power.t.test(n = df$n[i], sd = df$sd[i], delta = df$delta[i])$power
}
Data Science 2 | Folienübersicht

Sog. ‘Power Curves’

2

df %>% ggplot(aes(n, power, group = sd, color = factor(sd))) +
  geom_line() + geom_hline(yintercept = 0.8) + facet_grid(~delta, labeller = label_both)

  • Bei hoher Effektgröße und geringer Variabilität wird die nicht-lineare Beziehung zwischen der Stichprobengröße und der Power sichtbar.
    • → Kleine Erhöhungen des Stichprobenumfangs führen hier bereits zu großen Erhöhungen der Teststärke.
  • Die Kurve macht schnell deutlich, dass ein N von 40 nicht ausreicht um eine Power von 0.8 zu erzielen, wenn die Standardabweichung >750 ist und die Differenz bei nur 500km liegt.
Data Science 2 | Folienübersicht

Simulations-basierte Poweranalyse

1

  • Die Effektgröße und die Variabilität sind 2 Parameter in der Poweranalyse, die meist recht ungewiss sind.
  • Daher sind die berechneten Stichprobengrößen und Powerschätzungen mit hohen Unsicherheiten belegt.
  • Anstelle zu fragen “Welchen Stichprobenumfang brauche ich um 80% Teststärke zu haben?” sollte besser gefragt werden “Welchen Stichprobenumfang brauche ich um 90% sicher zu sein, dass das Experiment eine Power von 80% hat?”

Simulation für das Zugverhalten-Beispiel 

set.seed(123)
# Erzeugung von normalverteilten mittleren Differenzen und SDs
delta <- rnorm(5000, 1200, 100)
sigma <- rnorm(5000, 951, 100)
df <- data.frame(delta = delta, sd = sigma)

# Berechnung von N für jede Kombination an ES und SD
n <- apply(df, 1, function(x) power.t.test(delta = x[1], sd = x[2], power = 0.8)$n)
Data Science 2 | Folienübersicht

Simulations-basierte Poweranalyse

2

Wie hoch sollte N sein, um eine 90%ige Wahrscheinlichkeit zu haben, dass die Power 80% beträgt?

quantile(n, 0.9)
     90% 
14.91489

→ Unsere ursprüngliche Poweranalyse empfahl eine N von 11.

p
Data Science 2 | Folienübersicht

Berechnung des Stichprobenumfangs

ANOVA

power.anova.test(groups = NULL, n = NULL,
                 between.var = NULL, within.var = NULL,
                 sig.level = 0.05, power = NULL)

Input-Paramater

  • Anzahl der Gruppen
  • Stichprobenumfang pro Gruppe (Annahme balanciertes Design)
  • Varianz zwischen den Gruppen
  • Varianz innerhalb der Gruppen
  • Signifikanzniveau
  • Teststärke
Data Science 2 | Folienübersicht

Berechnung des Stichprobenumfangs

ANOVA - Beispiel

Gewichtsunterschiede von Atlantischem Lachs, der in 4 verschiedenen Typen von Netzkäfigen gezüchtet wurde (N = 24).

df <- salmon %>%  group_by(cages) %>% 
  summarise(mean = mean(weight), 
    sd = sd(weight))
power.anova.test(groups = 4, 
  between.var = var(df$mean), 
  within.var = mean(df$sd)^2, 
  power = 0.8)
     Balanced one-way analysis of variance power calculation 

         groups = 4
              n = 2.345194
    between.var = 0.7663649
     within.var = 0.2335063
      sig.level = 0.05
          power = 0.8

NOTE: n is number in each group
Data Science 2 | Folienübersicht

Power-Kurve

ANOVA - Beispiel

n <- 2:10
p <- power.anova.test(groups = 4, n = n, 
  between.var = var(df$mean), 
  within.var = mean(df$sd)^2)
p
     Balanced one-way analysis of variance power calculation 

         groups = 4
              n = 2, 3, 4, 5, 6, 7, 8, 9, 10
    between.var = 0.7663649
     within.var = 0.2335063
      sig.level = 0.05
          power = 0.6193382, 0.9544039, 0.9968026, 0.9998332, 0.9999929, 0.9999997, 1.0000000, 1.0000000, 1.0000000

NOTE: n is number in each group
plot(n,p$power)

Data Science 2 | Folienübersicht

Berücksichtigungen

  • Die Poweranalyse und Stichprobenkalkulation basieren auf den gleichen Annahmen wie die zugrunde liegenden Tests:
    • Für eine zuverlässige Poweranalyse und Stichprobenkalkulation beim t-Test muss die Normalitätsannahme für jede Gruppe erfüllt sein.
    • Bei der Poweranalyse für eine ANOVA müssen alle Gruppen dazu die gleiche gemeinsame Varianz haben, damit dem Ergebnis vertraut werden kann.
  • Da die Effektgröße und Variabilität oft unsicher sind, nutze die Simulations-basierte Analyse.
Data Science 2 | Folienübersicht

Hilfreiches Paket: pwr

  • Funktionen für weitere Tests und Design gibt es in dem R Paket pwr von Stéphane Champely, welches sich an Cohen (1988) orientiert (siehe auch die Vignette):
Funktion Poweranalyse für
pwr.2p.test() 2 Proportionen (gleiches N)
pwr.2p2n.test() 2 Proportionen (ungleiches N)
pwr.anova.test() Balanced one-way ANOVA (gleiches N)
pwr.chisq.test() Chi-Quadrat-Test
pwr.f2.test() Generalisierte Lineare Modelle (GLM)
pwr.p.test() Proportion (1 Stichprobe)
pwr.r.test() Korrelation
pwr.t.test() t-Tests (1-2-Stichproben, unabhängig, gepaart)
pwr.t2n.test() t-Test (2-Stichproben, ungleiches N)
Data Science 2 | Folienübersicht

Weitere hilfreiche Pakete

  • powerAnalysis: Funktionen für einfache Tests (ähnlich wie pwr)
  • webpower: Sammlung von Funktionen für simple Tests bis komplexe lineare Modelle, SEMs, etc.
  • Superpower: Poweranalysen für komplexere ANOVA Designs
  • simglm: Spezialpaket für Poweranalysen mittels Simulationen
  • SIMR, SIMR, SIMR: Poweranalyse für Generalized Linear Mixed Effects Models
Data Science 2 | Folienübersicht

Ãœbungsaufgabe

Data Science 2 | Folienübersicht

Übungen aus…

Kapitel 11 - Poweranalyse und Bestimmung des Stichprobenumfangs


  • R Notebook-Skripte
    • DS2_11_Ãœbungen.Rmd
    • DS2_11_Ãœbungen_Lösung.Rmd

Data Science 2 | Folienübersicht

Abschlussquiz

Data Science 2 | Folienübersicht

Fragen?

Data Science 2 | Folienübersicht